function [ T ] = enhanceGrid3D( T )
% function [ T ] = enhanceGrid3D( T )
%
%Input:
%     T : Extended tetrahedral data structure
%     
% Output:
%     T : Enhanced tetrahedral data structure with T.volume, T.normal,
%           T.area
%
% Last modified: Apr 10, 2015

nelts = size(T.elements,2);

T.volume=(1/6)*...
 ((T.coordinates(1,T.elements(2,:))-T.coordinates(1,T.elements(1,:))).*...
     (T.coordinates(2,T.elements(3,:))-T.coordinates(2,T.elements(1,:))).*...
     (T.coordinates(3,T.elements(4,:))-T.coordinates(3,T.elements(1,:)))-...
  (T.coordinates(1,T.elements(2,:))-T.coordinates(1,T.elements(1,:))).*...
     (T.coordinates(3,T.elements(3,:))-T.coordinates(3,T.elements(1,:))).*...
     (T.coordinates(2,T.elements(4,:))-T.coordinates(2,T.elements(1,:)))-...
  (T.coordinates(2,T.elements(2,:))-T.coordinates(2,T.elements(1,:))).*...
     (T.coordinates(1,T.elements(3,:))-T.coordinates(1,T.elements(1,:))).*...
     (T.coordinates(3,T.elements(4,:))-T.coordinates(3,T.elements(1,:)))+...
  (T.coordinates(2,T.elements(2,:))-T.coordinates(2,T.elements(1,:))).*...
     (T.coordinates(3,T.elements(3,:))-T.coordinates(3,T.elements(1,:))).*...
     (T.coordinates(1,T.elements(4,:))-T.coordinates(1,T.elements(1,:)))+...
  (T.coordinates(3,T.elements(2,:))-T.coordinates(3,T.elements(1,:))).*...
     (T.coordinates(1,T.elements(3,:))-T.coordinates(1,T.elements(1,:))).*...
     (T.coordinates(2,T.elements(4,:))-T.coordinates(2,T.elements(1,:)))-...
  (T.coordinates(3,T.elements(2,:))-T.coordinates(3,T.elements(1,:))).*...
     (T.coordinates(2,T.elements(3,:))-T.coordinates(2,T.elements(1,:))).*...
     (T.coordinates(1,T.elements(4,:))-T.coordinates(1,T.elements(1,:))));
 
T.area=(1/2)*sqrt(...
 ((T.coordinates(2,T.faces(2,:))-T.coordinates(2,T.faces(1,:))).*...
  (T.coordinates(3,T.faces(3,:))-T.coordinates(3,T.faces(1,:)))-...
  (T.coordinates(3,T.faces(2,:))-T.coordinates(3,T.faces(1,:))).*...
  (T.coordinates(2,T.faces(3,:))-T.coordinates(2,T.faces(1,:)))).^2+...
 ((T.coordinates(1,T.faces(2,:))-T.coordinates(1,T.faces(1,:))).*...
  (T.coordinates(3,T.faces(3,:))-T.coordinates(3,T.faces(1,:)))-...
  (T.coordinates(3,T.faces(2,:))-T.coordinates(3,T.faces(1,:))).*...
  (T.coordinates(1,T.faces(3,:))-T.coordinates(1,T.faces(1,:)))).^2+...
 ((T.coordinates(1,T.faces(2,:))-T.coordinates(1,T.faces(1,:))).*...
  (T.coordinates(2,T.faces(3,:))-T.coordinates(2,T.faces(1,:)))-...
  (T.coordinates(2,T.faces(2,:))-T.coordinates(2,T.faces(1,:))).*...
  (T.coordinates(1,T.faces(3,:))-T.coordinates(1,T.faces(1,:)))).^2);

% Normals to the faces
XX = reshape(T.coordinates(1,T.faces(1:3,:)),size(T.faces(1:3,:)));
YY = reshape(T.coordinates(2,T.faces(1:3,:)),size(T.faces(1:3,:)));
ZZ = reshape(T.coordinates(3,T.faces(1:3,:)),size(T.faces(1:3,:)));

VX1 = XX(2,:)-XX(1,:);
VY1 = YY(2,:)-YY(1,:);
VZ1 = ZZ(2,:)-ZZ(1,:);

VX2 = XX(3,:)-XX(1,:);
VY2 = YY(3,:)-YY(1,:);
VZ2 = ZZ(3,:)-ZZ(1,:);

V1 = [VX1;VY1;VZ1];
V2 = [VX2;VY2;VZ2];

T.normals = (1/2)*cross(V1',V2')';

return
